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Abstract 

In this paper we propose a new method for the study and visualization of 
dynamic processes in magnetic nanostructures, and for the accurate calcu- 
lation of rates for such processes. The method is illustrated for the case 
of switching of a grain of an exchange-coupled recording medium, which 
switches through domain wall nucleation and motion, but is generalizable to 
other rate processes such as vortex formation and annihilation. The method 
involves calculating the most probable (lowest energy) switching path and 
projecting the motion onto that path. The motion is conveniently visualized 
in a two-dimensional projection parameterized by the dipole and quadrupole 
moments of the grain. The motion along that path can then be described 
by a Langevin equation, and its rate can be computed by the classic method 
of Kramers (1940). The rate can be evaluated numerically, or in an analytic 
approximation - interestingly, the analytic result for domain-wall switching 
is very similar to that obtained by Brown in 1963 for coherent switching, 
except for a factor proportional to the domain wall volume. Thus in addi- 
tion to its lower coercivity, an exchange-coupled medium has the additional 
advantage (over a uniform medium) of greater thermal stability, for a fixed 
energy barrier. 
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1. Introduction 



Calculating switching rates of small magnetic elements is a very impor- 
tant problem in magnetic information storage. In addition to conventional 
field-switched media, new ideas for storing information in vortex chirality 
and polarity [l[ and spin torque switched elements j^j require calculation of 
switching rates, both with external fields (modeling the writing process) and 
without (modeling stability). 

Several approaches have been explored for this purpose, all ultimately 
involving the solution of a Fokker-Planck equation for the probability dis- 
tribution of the system. The problem was first approached by Brown [3[ in 
1963, who computed the switching rate for a coherent single-domain uni- 
axially symmetric particle. In this case the system is described by a single 
variable, the angle 9 of the magnetization from the easy axis, or alternatively 
the easy- axis component m x of the magnetic dipole moment. Once the prob- 
lem has been reduced to a one-dimensional one described by a stochastic 
Langevin-like equation of the form 



where the variance of the random term is a known function of m x , the rate can 
be determined by the well-established method of Kramers(1940) jij. Brown's 
approach can be generalized somewhat but is limited to coherently-switching 



If we go beyond coherent systems to a system with a spatially-varying 
magnetization described by a micromagnetic model with N cells, the config- 
uration space of the system is 2iV-dimensional and very hard to visualize. A 
general method for rate calculation, involving the expansion of the energy 
function in a power series near the saddle point in 2N-dimensional config- 
urational space (a Hessian matrix) was developed by Langer et al jsf and 
has been applied to magnetic systems However, this approach is com- 
putationally demanding (because it involves 2N x 2N matrices), somewhat 
formal (giving little heuristic insight into what physical factors control the 
rate), and because it involves only information near the saddle point, cannot 
be exact in the low temperature limit in the sense that the present approach 
appears to be. 

The approach we will take in the present paper is to map our 2N-dimensional 
system into a ID system governed by a Langevin equation, and compute the 
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switching rate of that ID system by the method of Kramers. We will develop 
the method in the context of an exchange-coupled or anisotropy-gradient el- 
ement, which is known [8] to switch incoherently by a sort of domain wall 
motion, from the initial state shown in Fig. [1] to a final reversed state (see 
Fig. [2] below). We use an anisotropy variation K(x) oc x 2 , which can be 
shown to give optimal switching in the limit of a long particle jjj. We also 
ignore magnetostatic fields in this calculation. Although it is very difficult 
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Figure 1: The 10 rim long anisotropy-graded grain considered in this paper. The anisotropy 
field Hk(x) oc x 2 , ranging from at the left (x = 0) to 9550 kA/m at the right (very 
high, but based on numbers for FePt). We use exchange constant A = 10~ n j/m and 
M s = 1520 kA/m. 

to make direct comparisons with conventional stochastic simulation (because 
simulation is very difficult for these very slow rates), it seems likely that this 
method is exact in the limit of low temperatures (large -Ebamer/ &bT), because 
the deviations from the minimum-energy path vanish in this limit. In this 
paper we will not deal with the question of global searches for saddle points 
and switching paths for arbitrary systems - we are using a relatively simple 
system (a graded-anisotropy grain) for which the switching mechanism (DW 
motion) is already known. 

2. Mapping to ID Langevin motion 

The switching involves motion along or near a one-dimensional path in 
2iV-dimensional space, the switching path, parameterized by a single real 
number, the switching coordinate. The simplest choice for this switching 
coordinate is of course m x , the moment along the easy (long) axis of the 
element, as in the original Brown theory. In this original theory, m x deter- 
mines the energy of the system, since the energy is independent of azimuthal 
angle. We could hope that even in an incoherently switching system, this 
single "coarse-grained" variable still essentially determines the state of the 
system, in the sense that if we specify m x , fluctuations in the other variables 
(except the irrelevant azimuthal angle) are small, because they are expen- 
sive in energy. That is, m x determines values for the other 2N — 2 variables 
(obtained by minimizing the energy at fixed m x ) from which fluctuations are 



3 



small, because the energy function rises rapidly when these other variables 
deviate from their minimum-energy values. This picture turns out to be true 
near the initial state, as a reversed domain forms, and as the domain wall 
moves. However, when this domain wall approaches its highest energy (when 
it is in the hard end of the system, see inset labeled "saddle point" in Fig. 
[2] below) m x no longer specifies the state of the system uniquely. There are 
many states with the same m x , and about the same energy. One way to see 
this is to observe that in this configuration, we can increase m x by moving 
the domain wall to the right (tilting the magnetizations near the domain 
wall counterclockwise), or by shortening the tail of the domain wall (tilting 
the magnetizations in the soft end counterclockwise). Near the saddle point, 
neither of these costs much energy. Doing these two things at the same time 
leaves 



unchanged. Thus a single coarse-grained variable is not sufficient to describe 
the system - we need to identify a second variable and hope that the state 
is determined by the two variables together. Noting that the two states of 
increased m x mentioned above have the increase in different places (in the 
hard end and the soft end respectively, i. e. at positive and negative x, if we 
take our origin at the center of the grain), we see that they will have different 
values of the weighted average 



Note that this is just the quadrupole moment. Thus we can hope to de- 
scribe the system by the two coarse-grained variables m x and q xx , in the 
sense that if we fix these two variables and minimize the energy with respect 
to the remaining 2N — 2 variables, the fluctuations of these remaining vari- 
ables (except, again, the azimuthal angle) are small. We can think of these 
remaining variables as "short wavelength spin waves". Denoting this min- 
imum value by E(m x ,q xx ), we can plot the contours of constant energy as 
functions of m x and q xx (Fig. [2]) In the context of this contour plot, we can 
understand more clearly the failure of m x to completely describe the system. 
If it did, the energy would rise rapidly as we move away from the minimum 
in any of the other 2N — 1 directions, including that described by q xx . That 
is, the energy at fixed m x would have a sharp minimum as a function of q xx . 
We can see in Fig. H]that this is true for the first part of the switching (at the 
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Figure 2: Contours of the energy function E(m x ,q xx ), showing also an approximate 
steepest-descent switching path. Insets show the actual minimum-energy configuration 
at various points in the m x -q xx plane. The boundaries of the colored region (maximum 
or minimum q xx ) are known analytically - they represent configurations with infinitely 
sharp domain walls (Fig. [3]). Contours are produced from data on a grid, by the Origin 
TM graphing p ac k a g e They cannot be accurate near the boundaries because the energy 
diverges in that limit, although they are accurate near the switching path. 
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Figure 3: System with infinitely sharp domain wall, represented by the bottom point in 
Fig. [3J with maximum \q xx \. 
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left) but the minima get shallower as we approach the saddle point, finally 
disappearing entirely at the saddle point. 




Figure 4: Same energy surface as in Fig. [2 with insets showing cross-sections through it 
at constant m x . Note that at the beginning of the switching process (on the left) there is a 
deep energy minimum, meaning that fluctuations about this configuration are small. Near 
the saddle point, this minimum becomes very shallow and disappears, indicating large 
fluctuations - specifying such a value of m x does not determine the state of the system. 
(When there are two minima, we are interested in the left-hand one - the other represents 
a reversed switching path at the top of the figure, related by M x <-> —M x ) 



So far, we have a 2D projection of the motion that allows us to visualize 
it, but we have not defined a unique switching path - we have seen that 
minimizing the energy at fixed m x is not a suitable method for doing so. 
Another way to define a switching path is by steepest descent in energy. 
At the saddle point, where 2N — 1 of the 2N principal curvatures of the 
energy function are positive, there is a unique direction in which the curvature 
is downward. Moving in this direction, and continuing along the energy 
gradient, defines a unique path which eventually reaches the minimum-energy 
initial saturated configuration (Fig. [1]). This appears to be the same path 



that is obtained by "nudged elastic band" methods [10|, lUl], but we have 
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computed it by following the energy gradient, which is equivalent to the use 
of the Landau-Lifshitz-Gilbert dynamical equation for the magnetization M 
of a single macrospin: 



(1 + a 2 ) 



dM 



(4) 



dt 



x H - — M x (M x H 



in the limit of infinite Landau-Lifshitz damping factor a, so it becomes 



We must let the gyromagnetic ratio 7 —¥ 00 as well, such that the ratio j/a 
remains finite. Physically, this eliminates the precession term and causes the 
system to move along the energy gradient in 2N-dimensional configuration 
space. 

Although the system cannot be described at all points of the switching 
process by specifying the single variable m x and allowing the other variables 
to relax to their minimum-energy values, this does not mean that the switch- 
ing process is not in some sense one-dimensional. The energy still goes up 
if we move in the direction normal to the steepest descent (SD) path (along 
the SD path is the only direction in which it goes down). If we define a 
new "switching coordinate" in the 2N dimensional space by by projecting 
each point perpendicularly onto the SD path, then fluctuations from these 
projected points are indeed small (the energy goes up rapidly perpendicular 
to this path). The system can still be described by a single variable, it just 
can't be m x in the sense that we described above. We can parameterize the 
path in any way we want - for example, we could use the path length along 
the SD path as a switching coordinate. However, this would be numerically 
complicated, and we could equally well use any monotonic function of this 
path length. Paradoxically, m x is such a monotonic function, though we just 
argued that it cannot be used to uniquely or almost-uniquely describe the 
system! The resolution of this paradox is that we cannot determine the other 
variables by minimizing the energy in a hyperplane of constant m x (a verti- 
cal plane in Fig. [2]) - there is no sharp minimum. However, if we minimize 
the energy in a hyperplane perpendicular to the SD path, there is a sharp 
minimum that determines all the other variables. 

This leads us to a prescription for mapping an arbitrary system trajectory 
onto a random walk along the SD path. At each point of the SD path, we 
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calculate the energy gradient (which is the tangent to the SD path) and draw 
the 2N — 1 dimensional hyperplane normal to this direction. If the trajectory 
is close to the SD path, it will be on one and only one of these hyperplanes, 
and this determines its switching coordinate. [If it is far away, farther than 
the radius of curvature of the SD path, it could be on two such hyperplanes 
- but this point would have very high energy, so we will assume that this is 
a rare event and worry about it later.] 

To describe the statistics of switching, we now must write an equation 
of motion for the switching coordinate, which should take the form of a 
Langevin equation. Since we are parameterizing the SD path by m x , this is 
an equation for dm x /dt. To avoid problems related to the Ito-Stratonovich 



ambiguity [12l. |13| . we find it better to discuss the change Am x during a 
finite time increment At of a micromagnetic simulation. We must be careful 
how we define the change Am x from the actual change in configuration (the 
magnetization change AM r at cell r). We cannot just substitute this change 
into Eq. [2] to get Am x , but we must first project the new configuration back 
to the SD path, and only then use Eq. [2J Being similarly careful about 
defining the diffusion term (variance of Am x ), we end up with a Langevin 
equation of the form of Eq. [TJ Since the dependence on the external field H x 
is linear, we may put it into the form 

Am x = fi(m x )[H x - H pin (m x )]At + (Am x ) random (6) 

This expression is exactly derived from the Landau-Lifshitz-Gilbert equation 
- however, it is important to realize that it is exact only on the ID SD 
path. To the extent that the real stochastic system deviates from this path, 
it is only approximate. The accuracy could be checked by examining the 
actual switching trajectories in a stochastic simulation, although this would 
be possible only for unrealistically low -Ebamer /&b T (in a real device, this 
should be at least 40). As T — )• 0, we expect the actual trajectories to become 
arbitrarily close to the SD path. We have defined the "field mobility" fi(m x ) 
as the coefficient of the external field in the function f(m x ) in Eq. [TJ it 
is a sum over the cells, whose most important characteristic is that it has 
contributions only from the domain wall, i. e., the contribution depends 
on the hard-axis component of the magnetization. In fact, an approximate 
expression for the field mobility is 

/i = Vwall (7) 

a 
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where V wa u is the volume of the domain wall (a small fraction of the total 
volume Vtotai)- We have also defined the pinning field (the field required to 
keep the domain wall stationary) 

H pm {m x ) ee (8) 
The random term is defined by its variance 

[(^random 2 ) = 2D(m x )At (9) 

which is a function of position m x along the path - D(m x ) can be calculated 
from the magnetization configuration {M r }. It is related to the mobility by 
an Einstein relation, 

kTjj(m x ) = n D(m x ) (10) 

This ID motion is slightly more general than Langevin's original formu- 
lation, in that the mobility and diffusivity (which are constant in a textbook 
Langevin equation) are functions of position in our Langevin equation. 

We use the method of Kramers[4| to compute the rate from our ID 
Langevin equation. This involves beginning with the Fokker-Planck equa- 
tion for the probability distribution of solutions to Eq. O We assume the 
distribution near the barrier is in a steady state, in which the probability 
current 

J = p/j(H x - H pin (m x )) - D^- (11) 

dm x 

is small and uniform. Following Kramers, we write the probability density 
function C(m x ) times its equilibrium value: 

p(m x ) = C(m x )e~ [E{mx) -^ H]/kBT (12) 

so that C(m x ) ~ in the switched well (m x +m s ), and is constant in 
the initial well. Inserting this into Eq. [TTJ letting H x = 0, and using the 
uniformity of the current gives an equation for this constant, 

C(-m s ) = r —L-e +E ^ k ° T dm x (13) 
J-m s D (m x ) 

The switching rate is then r = J/(well population), where in the populated 
part of the initial well, C(m x ) ~ C(—m s ), so 

J J 



HZ! P{m x )dm x C(-m s )f+™ s e~ E( - m ^ k » T dm, 
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J+^ s s D(m x )- 1 e+ E ( m ^/ k B T dm x f+^e- E ( m -)/ k B T dm x v ; 

The first integral in the denominator is dominated by the vicinity of the 
barrier, where E(m x ) pa -Ebamer- We can bring out the Arrhenius factor to 
give 

where the prefactor ("attempt frequency") is 

r ° = /+™; D(m x )- l e+y E ^)- E ^y k BTdm x /+™; e- E ( m *)l k BTdm x ^ 

Using Eq. [161 we can numerically calculate the integrals appearing in 
the rate prefactor from our calculations of the states along the SD (steepest- 
descent) path. However, it is also interesting to look at the effects of some 
simple approximations to these integrals, because it sheds light on the con- 
nection with the Brown result for coherent switching. Brown's 1963 result 
for the prefactor of a uniaxial particle with anisotropy field Hk was [3] 

r °=w^r (i7> 

We can put our result in a similar form - first note that the first integral in the 
denominator of Eq. [T3]is very sharply peaked at the energy barrier, especially 
at low temperatures. We can approximate the energy by a quadratic function 
near its peak 

F(m)-F (\ K-^^H /io\ 

L{m x ) - ^barrier I 1 I (18) 

where B is a dimensionless barrier width factor, of order 1. Also at reasonably 
low temperatures, the second integral in the denominator of Eq. [TH depends 
only on the initial slope of E(m x ) near the initial state, which is the nucleation 
field H nnc . 

These approximations allow us to evaluate both integrals analytically, 
giving us 

r = a 7 H nuc (—±-) ^ (19) 



.7rk B Tj B 

Surprisingly, this is identical to Brown's result for the coherent case, except 
that the nucleation field is playing the role of the anisotropy field, and there 
are factors of the domain wall fraction W and the barrier width parameter 
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B, which should be of order 1. The proportionality to the domain wall width 
is easy to understand physically, because magnetization fluctuations outside 
the domain wall cannot drive domain wall motion. 

To check the importance of continuous anisotropy grading (as opposed to 
having layers of uniform but different anisotropy) we considered a series of 
systems obtained from a continuously graded system by constructing layers 
and assigning each a uniform anisotropy equal to the average over its volume 
in the continuously graded system. As the number of layers increases, we thus 
expect the behavior to get closer to that of the continuous system. However, 
with 1 layer (a uniform system) we expect coherent rotation, then with 2 
or more layers the switching mechanism involves a domain wall. From Eq. 
dHJ we expect a large decrease in rate prefactor when a domain wall appears 
(W « 1), and this is confirmed by the numerical results (Fig. [5]). There is 
suprisingly little change if we increase the number of layers beyond 2. The 
advantage of using more than 2 layers lies instead in the improved coercivity. 
We have also computed the coercivity as a function of number of layers, 
and find that the coercivity increases as we decrease the number of layers, 
because this increases the anisotropy jump at the interfaces and therefore 
the strength of the domain wall pinning at these interfaces. 

3. Conclusion 

In this paper we have outlined a new method for the calculation of switch- 
ing rates of incoherently switching micromagnetic systems, and applied it to 
a simple graded-anisotropy system. The result can be put into a form that 
is surprisingly similar to the form derived in 1963 by Brown for coherent 
switching, except for factors that can be physically interpreted as effects of 
incoherence, such as the domain wall width. The predicted rates can be 
tested by comparison with direct simulation, and such tests are planned. 
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Figure 5: Thermal switching rate prefactor as function of coarseness of discretization 
(number of layers) . The data with black squares are obtained by parabolic fitting of the 
energy landscape (Eq. I18[) . The red circles are obtained by direct numerical integral of the 
energy landscape (Eq. fT6|) . The switching rate prefactor of the Stoner-Wohlfarth particle 
(one layer) is obtained by analytical calculation of the energy landscape assuming that the 
energy barrier is the same as a 20-cell anisotropy-graded system. 
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